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Accurate  and  numerically  efficient  modeling  of  low  to  moderate  Reynolds 
number  nozzle  flow  expansions  to  vacuum  can  be  difficult  due  to  the  pres¬ 
ence  of  multiple  flow  length  scales.  Such  simulations  are  important  for  the 
prediction  of  propulsive  thrust  as  well  as  spacecraft  contamination,  both  of 
which  can  be  difficult  to  measure  in  ground  based  facilities.  To  that  end, 
conical  nozzle  flows  were  studied  for  Reynolds  nnmbers  of  1,230  and  12,300 
using  the  direct  simnlation  Monte  Carlo  method  (DSMC),  Navier-Stokes 
with  velocity  slip  and  temperature  jump  boundary  conditions,  and  statisti¬ 
cal  and  deterministic  approaches  to  the  solntion  of  the  BGK  and  ES-BGK 
equations.  The  deterministic  and  statistical  solutions  of  the  BGK  equation 
were  found  to  be  in  good  agreement  with  the  benchmark  DSMC  results. 
Statistical  BGK  and  ES-BGK  methods  were  also  found  to  be  more  efficient 
methods  than  DSMC  in  the  continunm  and  near-continunm  regime,  and 
more  accurate  than  the  Navier-Stokes  eqnations  in  the  portions  of  the  flow 
with  rarefaction,  such  as  the  boundary  layer  and  the  flow  around  the  nozzle 
lip. 
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*Nomenclature 

Kronecker  delta 

Dynamic  viscosity 
Reference  dynamic  viscosity 
Characteristic  relaxation  frequency 
Density 
Solid  angle 

Viscosity  temperature  index 
Differential  cross-section  of  the  binary  collision 
Nozzle  throat  diameter 
Finite  volume 

Single  particle  velocity  distribution  function 
Maxwellian  distribution  function 
Ellipsoidal-Statistical  distribution  function 
External  force  per  unit  mass 
Boltzmann  constant 
Mass  of  a  molecule 
Number  of  particles  in  a  cell 

Number  of  particles  selected  for  velocity  resampling 
Number  density 
Prandtl  number 
Internal  pressure 
Radial  vector 
Random  numbers 
Reynolds  number 
Modifying  tensor 
T  :  Local  temperature  in  a  cell 

To  :  Stagnation  temperature 

Tref  ■  Reference  temperature 

t  :  Time 

U*  :  Sonic  velocity  at  the  nozzle  throat 

V  :  Velocity  vector 

Vr  '■  Relative  velocity  of  the  colliding  particles 

Vj  :  Velocity  vector  in  Cartesian  coordinates 
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I.  Introduction 


Nozzle  flows  at  low  and  moderate  Reynolds  numbers  are  characterized  by  multiple  flow 
length  scales,  which  signihcantly  complicates  their  accurate  and  numerically  efficient  mod¬ 
eling.  Multiparameteric  analysis  of  such  flows  is  also  problematic  due  to  a  large  number  of 
cases  that  need  to  be  examined  for  the  optimization  of  the  nozzle  based  device  performance. 
Experimental  studies  under  such  conditions  are  rare,  expensive  (particularly  for  microscale, 
MEMS  flows),  and  may  not  provide  the  necessary  accuracy  in  the  measurement  of  the  prin¬ 
cipal  nozzle  characteristics  such  as  thrust,  flow  rate,  and  specihc  impulse.^  For  both  micro 
and  meso-sized  nozzles  operating  in  the  space  near- vacuum  environment,  the  interaction 
with  critical  spacecraft  surfaces  at  high  altitudes  needs  to  be  analyzed.  Back-flow  produced 
by  such  devices  plays  a  major  role  in  the  contamination  to  sensitive  electronic  devices  such 
as  optical  instruments  and  solar  panels,  which  in  turn  may  adversely  affect  the  life  span  of 
spacecraft. The  back-flow  formation  is  very  sensitive  to  the  conditions  at  the  nozzle  lip, 
and  is  difficult  to  study  in  ground  based  facilities,  particularly  for  micro-nozzle  flows. ^  For 
the  International  Space  Station,  the  control  thruster  exhaust  products  turn  out  to  be  the 
most  important  agents  of  contamination.  The  contamination  of  such  electronic  devices  on  a 
spacecraft  can  become  mission  critical  by  signihcantly  reducing  its  active  life  span.^  There¬ 
fore,  in  the  present  work,  we  carried  out  nozzle  supersonic  expansion  to  vacuum  studies  with 
the  goal  of  studying  the  core  how  as  well  as  how  reversal  rate  diherences  among  the  statistical 
and  deterministic  ES-BGK  and  Navier-Stokes  (NS)  methods.  The  baseline  direct  simulation 
Monte  Carlo  (DSMC)  method  is  utilized  as  the  basis  for  comparison  for  the  aforementioned 
gas  dynamic  approaches.  The  numerical  accuracy  can  be  well  characterized  for  DSMC,  as 
will  be  discussed,  but  it  is  too  computationally  intensive  for  the  range  of  conditions  necessary 
to  analyze  micro-nozzle  thrust  and  spacecraft  contamination.  Therefore,  in  the  comparison 
among  the  gas  dynamic  techniques,  computational  efficiency  will  also  be  addressed. 

The  development  of  accurate  numerical  tools  capable  of  simulating  micro-nozzle  hows 
is  important,  but  at  the  same  time,  challenging,  because  the  how  regime  changes  from 
continuum,  near  the  nozzle  throat,  to  transitional  at  the  nozzle  exit.  Kinetic  methods,  such 
as  DSMC,  and  continuum  techniques  based  on  the  solution  of  the  NS  equations,  encounter 
computational  constraints  and  physical  challenges  when  applied  to  these  hows.  The  principal 
problem  with  the  DSMC  method^  is  the  associated  computational  cost  when  high  density 
portions  of  the  how  have  to  be  accurately  modeled.  This  is  essentially  due  to  the  large 
number  of  collisions  in  the  high  density  regime,  in  addition  to  the  usual  constraints  of  cell 
size  and  time  step  being  of  the  order  of  mean  free  path  and  mean  collision  time,  respectively. 
On  the  other  hand,  conventional  continuum  CFD  techniques  are  inapplicable  in  how  regions 
of  high  gradients  and  strong  rarefaction  even  when  velocity  slip  and  temperature  jump 
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boundary  conditions  are  imposed  at  the  nozzle  surface.  These  factors  become  paramount 
at  the  exit  of  the  nozzle  for  supersonic  expansions  to  vacuum,  but,  can  even  be  seen  in 
supersonic  expansions  to  non- vacuum,  hnite  back  pressures.®  Attempts  to  correct  the  NS 
equations  through  the  use  of  the  Burnett  equations  for  micronozzle  flow  modeling  showed 
that  the  latter  are  in  better  agreement  with  DSMC,  but  only  up  to  Knudsen  numbers  of 
approximately  0.2.^ 

In  recent  years,  combined  but  uncoupled  NS/DSMC  approaches  have  been  often  used  to 
model  nozzle  flows  in  the  transitional  Knudsen  number  regime,®’®  where  the  NS  equations 
were  used  to  model  the  high  density  portions  of  the  flow  inside  the  nozzle,  and  the  DSMC 
method  was  used  to  predict  the  rarehed  plume-atmosphere  interaction.  A  starting  surface 
was  used  to  transfer  the  properties  obtained  in  the  continuum  portion  of  the  flow  to  the 
DSMC  solver.  The  flow  was  always  assumed  to  be  in  the  steady  state,  and  the  uncoupled 
approach  implied  that  the  more  rarehed,  DSMC  portions  of  the  how  have  negligible  ehect 
on  the  continuum  portions.  This  approach,  however,  becomes  problematic  when  the  implicit 
assumptions  are  violated,  for  example,  when  the  hand-oh  surface  is  in  a  subsonic  region,  or 
the  downstream  how  ahects  the  upstream  how  through  radiation,  or  the  how  geometry  is 
three-dimensional  and  complicated.  La  Toree  et  recently  used  DSMC,  NS,  and  coupled- 
hybrid  DSMC/NS  approaches  to  evaluate  micronozzle  performance  for  expansions  to  vacuum 
and  found  the  macroparameter  distributions  to  be  sensitive  to  the  location  of  the  interface 
between  the  two  methods. 

To  overcome  this  problem,  it  is  possible  to  use  a  two-way  coupled  hybrid  continuum- 
kinetic  approach  (see,  for  example.  Refs.  11-13),  where  the  interface  location  can  be  dynam¬ 
ically  varied  during  the  simulation  by  evaluating  continuum  breakdown  parameters.  How¬ 
ever,  an  accurate  determination  of  the  interface  boundary  demands  additional  computational 
ehort.  There  will  also  be  the  usual  hybridization  problems^"^  and,  moreover,  the  extension  to 
chemically  reacting  hows,  i.e.,  to  hows  consisting  of  species  with  signihcant  internal  energy, 
or  to  two-phase  models  may  not  be  straightforward.  The  dehciencies  in  the  implementation 
of  a  hybrid  NS-DSMC  approach  become  even  more  problematic  when  transient  nozzle  and 
plume  hows  need  to  be  analyzed.  The  main  difficulty  is  related  to  the  temporal  changes  in 
gas  properties:  most  importantly  gas  mean  free  path  and  the  Mach  number,  that  necessitate 
hexible,  transient  hand-oh  surfaces  and  interface  boundaries. 

Considering  the  aforesaid  difficulties  in  modeling  nozzle  hows  in  the  transition  regime, 
it  would  be  desirable  to  have  a  single  method  that  allows  an  accurate  and  efficient  one-step 
modeling  of  high  density  nozzle  and  low  density  plume  hows.  Recent  years  have  witnessed 
a  renewed  interest  and  signihcant  advances  in  the  solution  of  model  kinetic  equations  such 
as  BGK,^®  BGK/ES-BGK,^®  with  deterministic,  either  hnite  diherence  or  hnite  volume, 
approaches  typically  used  in  the  solution  procedure.  A  particle  approach  to  obtain  the 
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solution  to  the  BGK  equation  was  first  proposed  in  Ref.  17.  It  was  then  extended  to  model  the 
ES-BGK  equation  in  Ref.  16,  and  further  extended  to  include  rotational  degrees  of  freedom 
in  Ref.  18.  The  general  idea  behind  the  use  of  model  kinetic  equations  as  replacement  for 
DSMG  is  that  the  solution  of  the  former  will  be  more  efficient  than  DSMG  in  the  continuum 
and  near-continuum  regimes,  and  more  accurate  than  the  solution  of  the  NS  equations  in 
the  transitional  regions  of  the  flow. 

With  respect  to  nozzle  flows,  one  of  hrst  attempts  to  use  model  kinetic  equations  is  the 
work  of  Burt  et  where  a  statistical,  particle  approach  to  the  solution  of  the  ES-BGK 
equation  was  obtained  for  a  flow  expanding  through  a  conical  nozzle.  More  recently  their 
work^®  was  extended  to  combine  the  advantages  of  the  NS  and  DSMG  approaches.  In  the 
low  diffusion  particle  method  the  viscosity  is  modeled  by  computing  the  diffusive  terms  in 
the  compressible  NS  equations  in  the  equilibrium  regions  of  the  flow  whereas  the  DSMG  pro¬ 
cedure  is  used  in  the  nonequilibrium  regions.  Although  initially  successful,  the  rehnement 
of  the  method  to  realistic,  more  complex  flows  including  internal  degrees  of  freedom  and 
chemistry  still  remains.  In  our  earlier  work,^*^  we  developed  and  studied  a  statistical  tech¬ 
nique,  that  models  continuum  flows  using  a  collision  enforcement  procedure  to  guarantee  full 
relaxation  of  the  particle  thermal  velocities  to  a  state  of  local  equilibrium,  hence  the  name, 
“eDSMG.”  The  technique  was  applied  to  internal  flows  for  nozzles  as  well  as  channels  and 
compared  well  with  nozzle  data  and  previous  channel  calculations.  The  technique  was  shown 
to  be  able  to  solve  inviscid  flows  with  tangency  (or  specular)  wall  boundary  conditions  {i.e., 
Eulerian),  but  under  predicted  viscous  effects  in  the  boundary  layer  when  diffuse  gas- wall 
boundary  conditions  were  used. 

In  the  present  work,  we  continue  our  efforts  to  apply  statistical  methods  to  moderate 
and  high  Reynolds  number  nozzle  flows  by  making  use  of  the  BGK^^  and  ES-BGK  models. 

A  particular  goal  is  to  understand  the  benehts  of  the  particle  BGK/ES-BGK  approaches 
as  compared  to  the  conventional  DSMG  method.  For  BGK/ES-BGK  approaches  to  be 
competitive  with  DSMG,  it  must  be  understood  if  the  standard  DSMG  requirements  can  be 
relaxed.  If  so,  this  would  then  allow  BGK/ES-BGK  models  to  be  used  for  high  pressure 
cases,  where  the  DSMG  method  is  impractical  due  to  the  required  computational  resources. 
As  it  will  be  shown  in  Sec.  nm  the  use  of  NS  method  is  not  a  good  alternative  for  supersonic 
expanding  flows  to  vacuum.  We  compare  numerical  simulations  for  transitional  nozzle  flows 
expanding  into  a  vacuum  using  Eve  different  numerical  approaches  -  DSMG,  hnite  volume 
and  particle  solutions  of  the  BGK  and  ES-BGK  model  kinetic  equations,  solution  of  the  NS 
equations  and  an  equilibrium  DSMG  (eDSMG)  technique  for  two  test  cases  corresponding 
to  Reynolds  numbers  of  1,230  and  12,300.  Such  comparisons  of  the  model  kinetic  equation 
approaches  with  the  formally  derived  computational  schemes  such  as  DSMG  and  NS  are 
rarely  found  in  a  single  article  and  will  allow  us  to  assess  the  utility  of  the  particle  BGK  and 
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ES-BGK  approaches  as  stand  alone  procednres.  Moreover,  in  futnre  work  it  may  be  possible 
to  ntilize  the  best  possible  featnres  of  the  particle-statistical  BGK  and  DSMG  methods  in 
terms  of  numerical  efficiency  and  accuracy  to  construct  a  particle-particle  hybrid  method 
consisting  of  these  two  methods. 

The  remainder  of  the  paper  is  as  follows.  In  Sec.  in]we  outline  the  theory  and  our  method 
of  implementation  of  the  computational  approaches  that  are  used.  We  discuss  the  numerical 
parameters  selected  for  the  DSMG  calculations  since  this  baseline  method  will  be  used  to 
model  a  high  pressure  case,  usually  outside  the  typical  practical  region  for  most  DSMG 
calculations.  In  Sec.  nm  we  discuss  the  numerical  parameters  chosen  for  the  BGK/ES-BGK 
calculations  to  minimize  statistical  errors  so  that  small  differences  among  the  methods  can 
be  seen.  Macroparameters  such  as  gas  density,  temperature,  and  velocity  along  the  nozzle 
centerline  and  across  the  nozzle  exit  are  compared,  as  well  as  back  flow  regions  downstream  of 
the  nozzle  exit.  To  compare  the  numerical  efficiency  of  the  methods,  the  rationale  for  a  new 
set  of  numerical  parameters  is  discussed  and  employed.  In  addition,  boundary  layer  growth 
and  flow  averaged  quantities  such  as  thrust  and  back  flow  rates  are  compared.  Gonclusions 
are  presented  in  Sec.  |IVl 


II.  Computational  Methods 
II.A.  DSMG  Method 

The  simulation  of  rarehed  flows  is  typically  performed  using  the  DSMG  approach,  a  dis¬ 
crete  particle  simulation  method  that  approximates  a  numerical  solution  to  the  Boltzmann 
equation.  DSMG  achieves  computational  efficiency  by  decoupling,  in  time,  the  movement 
and  collisions  of  molecules.  Simulated  particles,  each  representing  a  large  number  of  real 
atoms/molecules,  are  moved,  tracked  through  the  computational  mesh,  and  are  indexed  and 
sorted  into  cells  according  to  their  locations.  During  each  time  step,  some  fraction  of  the  par¬ 
ticles  in  a  cell  collide  with  each  other,  and  probabilistic  techniques  are  used  to  evaluate  the 
outcome  of  individual  collisions.  The  colliding  particles  are  assigned  new  velocities  accord¬ 
ing  to  the  specihed  interaction  law,  so  that  the  energy  and  momentum  are  conserved,  and 
are  then  moved  through  the  mesh  at  the  assigned  velocities.  Boundary  conditions  are  mod¬ 
eled  through  collisions  of  the  computational  particles  with  the  surfaces  by  properly  choosing 
the  gas-surface  interaction  parameters.  Macroparameters  are  obtained  by  the  appropriate 
summation  over  all  the  particles  contained  in  the  cells. ^ 

The  DSMG-based  SMILE  (Statistical  Modeling  in  a  Low  Density  Environment)  compu¬ 
tational  solver  is  used  in  this  work.  The  solver  has  both  2D  and  axisymmetric  capabilities. 
The  important  features  that  are  used  in  the  present  work  include  the  axisymmetric  capabil¬ 
ity  with  radial  weights,  different  grids  for  collisions  and  macroparameters,  both  of  which  are 


Distribution  A:  Approved  for  public  release;  distribution  unlimited 


two-level  adaptable  Cartesian  grids,  and  parallel  implementation  with  efficient  load  balanc¬ 
ing  techniqnes.^^  The  physical  space  is  modeled  in  a  cylindrical  coordinate  system,  while  the 
velocity  space  is  3-D  Cartesian.  The  majorant  frequency  scheme  of  Ivanov  and  Rogasinsky^^ 
is  employed  for  modeling  gaseous  species  collisions  and  the  Variable  Hard  Sphere  (VHS)^ 
model  is  used  for  modeling  the  total  collision  cross-section.  Diffuse  reflection  with  full  ther¬ 
mal  accommodation  is  assumed  at  the  nozzle  wall.  For  both  Reynolds  numbers,  solutions 
independent  of  grid,  time  step,  and  number  of  particles  are  obtained.  For  Re=l,230  (Case 

I,  see  Table  [T},  there  is  virtually  no  difference  observed  between  solutions  obtained  for  0.47 
million  cells  with  5  million  simulated  molecules,  and  3  million  cells  with  30  million  simulated 
molecules.  For  Re=12,300  (Case  II,  see  Table  El),  this  is  true  for  numerical  parameters  up 
to  an  order  of  magnitude  larger.  Tables  ^  and  El  provide  a  summary  of  the  numerical  pa¬ 
rameters  used  for  the  two  cases  considered  in  this  work.  The  DSMC  and  BGK  simulations 
are  accurate  to  within  2%.  The  eDSMC  calculations  discussed  in  Sec.  ehi  were  computed 
using  0.85  million  cells  with  5  million  particles. 

II.  B.  Solution  of  the  NS  Equations 

The  commercial  code,  CFD-|--|-,^'^  has  been  used  in  this  work  to  solve  the  NS  equations. 
CFD-I--I-  is  a  flexible  computational  fluid  dynamics  software  suitable  for  the  solution  of  the 
steady  and  unsteady,  compressible  and  incompressible  NS  equations,  and  includes  a  multi¬ 
species  capability  for  perfect  and  reacting  gases.  In  this  work,  a  perfect-gas  compressible  NS 
solver  is  used  with  second  order  spatial  discretization  and  implicit  time  integration.  Second 
order  velocity  slip  and  temperature  jump  conditions  are  imposed  on  the  nozzle  wall.  A  su¬ 
personic  inflow  with  prescribed  parameters  is  applied  at  the  nozzle  throat,  and  backpressure 
of  1  Pa  is  imposed  at  the  outflow  boundaries.  A  symmetry  condition  is  defined  at  the  nozzle 
axis.  The  results  presented  below  (both  flow  fields  and  computational  requirements)  are 
obtained  for  a  multi-block  rectangular  grid  with  a  total  of  14,400  nodes.  The  computations 
are  also  conducted  for  four  times  smaller  and  four  times  larger  numbers  of  nodes,  and  found 
to  be  fully  grid-resolved  with  14,400  nodes. 


II. C.  Solution  of  the  BGK  Equation 

The  velocity  distribution  function  provides  a  full  description  of  a  gas  at  the  molecular  level. 
The  relationship  for  the  velocity  distribution  function  is  given  by  the  Boltzmann  equation: 
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where  /  is  the  single  particle  velocity  distribution  function,  n  is  the  number  density,  F  is  an 
external  force  per  unit  mass  that  applies  to  the  particles  (assumed  to  be  zero  for  the  present 
study).  The  term  on  the  right  hand  side  is  the  collision  term  and  is  given  by  the  following 
equation: 


/oo  nAir 

/  irfi  -  ffi)  VrcrdVtdv  (2) 

■OO  »/0 

where  /  and  /i  are  the  values  of  the  velocity  distribution  function  at  pre-collision  velocities 
of  the  two  colliding  particles  and  f*  and  are  the  corresponding  values  at  post-collision 
velocities;  a  is  the  differential  cross-section  of  the  binary  collision,  is  the  solid  angle  and 
Vr  is  the  relative  velocity  of  the  colliding  particles.  It  is  clear  that  the  collision  term  involves 
multiple  integrations  in  its  analytical  formulation,  and  is  therefore  difficult  to  compute. 
Hence  different  simplihed  models  have  been  introduced  to  model  the  inherently  complicated 
collision  term  of  the  Boltzmann  equation.  One  such  simplihed  model  is  from  Bhatnagar, 
Gross  and  Krook,^^  that  approximates  the  collision  term  as  follows: 


dt 


collision 


^n{fe  -  f) 


(3) 


where  n  is  the  number  density,  u  is  the  characteristic  relaxation  frequency  and  fe  is  the 
Maxwellian  distribution  function.  In  this  model,  the  nonlinear  collision  term  in  the  Boltz¬ 
mann  equation  is  approximated  by  a  Maxwellian  molecular  model  in  which  the  distribution 
function  /  proceeds  toward  local  equilibrium  at  a  velocity  independent  rate.  Thus  the  BGK 
method  provides  an  alternate  procedure  to  account  for  the  collisional  process  driving  a  how 
toward  equilibrium  without  modeling  individual  collisions.  The  idea  behind  this  simplihca- 
tion  is  that  the  intricate  details  of  the  two-body  interactions  are  not  important  in  reproducing 
most  of  the  experimentally  measured  macroscopic  quantities  such  as  temperature,  pressure 
or  how  velocity,  if  the  collision  rate  is  sufficiently  high.^® 

The  BGK  equation  reproduces  correct  moments  and  satishes  the  H-theorem  for  en¬ 
tropy  production.  In  Eq.  El  the  term  nufe  represents  the  collisions  replenishing  the  local 
Maxwellian  equilibrium  distribution  fe,  and  the  term  nuf  represents  the  collisions  depleting 
the  molecules  out  of  existing  distribution  There  are  arguments  about  the  importance 
and  the  specihcs  of  algorithms  to  preserve  energy  and  momentum  conservation,  which,  in 
the  opinion  of  the  authors  is  a  secondary  problem,  since,  on  average,  these  quantities  are 
preserved.  The  major  issue  seems  to  be  the  absence  of  the  formal  derivation  of  the  nu¬ 
merical  schemes  from  the  BGK  equation,  in  a  similar  manner  as  it  is  conducted  for  DSMG 
schemes. For  this  reason,  the  simulation  results  of  BGK  methods  presented  in  Sec.  Illllwill 
be  compared  with  the  DSMG  method.  The  BGK  equation  is  solved  by  both  statistical  and 
hnite  volume  methods  in  the  present  work. 
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II.C.l.  Statistical  Method  for  Solution  of  the  BGK  Equation 

Recently,  a  number  of  authors^®’ have  developed  particle  approaches  to  the  solution  of 
the  BGK  equations.  While  these  approaches  differ  in  details  and  the  arguments  regarding 
the  formal  connection  between  the  numerical  schemes  and  the  equations  continue,  the  basic 
ideas  of  the  numerical  schemes  remain  the  same.  Instead  of  selecting  collision  pairs  from 
all  of  the  simulated  particles,  a  fraction  of  particles  are  randomly  selected  during  each  time 
step  from  those  available  in  a  computational  cell,  and  are  assigned  new  velocities  according 
to  local  Maxwellian  (BGK)  or  ellipsoidal  (ES-BGK)  distribution  functions.  If  the  collision 
frequency  for  such  a  velocity  reassignment  is  properly  computed  and  the  local  values  of  the 
translational  temperature  in  the  cells  are  known,  then  the  procedure  certainly  mimics  the 
collision  term  on  the  right  hand  side  of  the  BGK  equation  given  by  Eq.  El 

Statistical  BGK  Method 

Since  the  only  difference  between  the  DSMG  method  and  a  particle  solution  of  the  BGK 
equation  lies  in  the  modeling  of  the  collision  term,  procedures  other  than  intermolecular 
collisions  in  the  BGK  particle  method  are  identical  to  those  in  DSMG  method.  The  details 
of  the  statistical  BGK  model  are  as  followsd®^^® 

The  characteristic  relaxation/collision  frequency  is  calculated  as 

u  =  FT-nk  (4) 

Kl^ref  J 

where  Pr  is  the  Prandtl  number  (1  for  the  BGK  equation),  k  is  the  Boltzmann  constant,  T 
is  the  local  translational  temperature  in  a  cell,  oo  is  the  viscosity  temperature  index  and  p^e/ 
is  the  gas  dynamic  viscosity  at  Tref-  The  collision  frequency  is  calculated  for  each  compu¬ 
tational  cell  at  each  time  step  based  on  the  local  translational  temperature  T  and  the  local 
number  density  n.  The  local  number  density  n  is  averaged  over  a  large  number  of  compu¬ 
tational  time  steps,  while  the  local  temperature  T  is  computed  based  on  the  instantaneous 
thermal  velocities  of  the  computational  particles  in  the  cell. 

The  number  of  particles  preselected  for  velocity  resampling,  i.e.,  particles  for  which  ve¬ 
locities  are  reassigned  according  to  the  local  Maxwellian  distribution,  is  calculated  as  follows: 

Nc  =  int(iV(l  —  exp(— z/At))  (5) 

where  N  is  the  number  of  particles  in  a  cell.  At  is  the  time  step,  and  int  operator  means 
the  nearest  smaller  integer.  To  compensate  for  the  systematic  error  that  such  an  operator 
produces,  one  more  particle  is  added  to  to  the  list  of  preselected  particles  with  the  probability: 
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Pc  =  N{1  -  exp(-i/At)  -  int(A(l  -  exp(-i/At))  (6) 

The  preselected  particles  receive  new  velocities,  ,  according  to  the  local  Maxwellian  dis¬ 
tribution.  Each  velocity  component  is  set  via 

vl  =  cos(27rPi) y/— ln(P2)  •  \/2kT/m  (7) 


Vy  =  sin(27rPi) ln(P2)  •  ^/2kT/m 


(8) 


vl  =  cos(27rP3)  y/—  ln(P4)  •  \/2kT/m  (9) 

where  Pi  through  P4  are  random  numbers  uniformly  distributed  between  0  and  1.  The  ve¬ 
locities  of  particles  that  have  not  been  preselected  remain  unchanged  in  the  current  time  step. 


Statistical  ES-BGK  Method 

The  BGK  model  has  an  inherent  property  that  Prandtl  number  is  unity  for  all  cases.  This 
may  be  a  source  of  error  for  cases  in  which  the  thermal  conductivity  plays  an  important  role. 
Holway^®  and  Cercignani^^  proposed  a  modification  to  the  BGK  equation  that  ensures  the 
correct  production  of  transport  coefficients  similar  to  the  Boltzmann  equation,  and  has  been 
recently  shown  to  satisfy  Boltzmann’s  H-theorem  for  entropy  production.^®  The  ES-BGK 
model  is  obtained  by  the  replacement  of  the  Maxwellian  distribution  with  a  local  anisotropic 
three-dimensional  Gaussian,  referred  to  as  the  Ellipsoidal-Statistical  (ES)  model: 


dt 


collision 


vn(  fa  -  f) 


(10) 


The  procedure  of  calculating  the  characteristic  relaxation  frequency  (with  the  physically 
correct  Prandtl  number),  selection  and  resampling  of  velocities  for  a  fraction  of  particles 
available  in  a  computational  cell  are  the  same  as  in  the  statistical  BGK  method.  However 
after  resampling  of  all  velocity  components,  ,  they  are  modified  to  in  order  to  conform 
to  the  ES  distribution; 

Vi=Sij-vj  (11) 

where  designates  the  modified  velocity  components  and  Sij  is  given  by  the  following 
equation: 
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where  the  symbol  <>a  represents  an  averaging  over  all  the  particles  in  a  cell,  dij  is  the 
Kronecker  delta,  and  m  is  the  mass  of  a  particle. 

II. D.  Finite  Volume  Method  for  Solution  of  the  BGK/ES-BGK  Equation 

A  hnite  volume  solver  SMOKE  developed  at  ERG,  Inc.,  has  been  used  to  deterministically 
solve  the  BGK  and  ES-BGK  equations.  SMOKE  is  a  parallel  code  based  on  numerical 
schemes  developed  by  L.  Mieussens.^®  Both  implicit  and  explicit  time  integration  schemes 
are  implemented.  The  time-accurate  explicit  schemes  are  constructed  to  strictly  satisfy 
mass,  momentum,  and  energy  conservation  requirements,  while  the  implicit  schemes  do  not 
have  this  advantage.  Nevertheless,  for  the  problem  under  consideration  it  was  found  that 
the  implicit  schemes  provide  adequate  accuracy,  and  subsequent  application  of  the  explicit 
schemes  essentially  does  not  visibly  change  the  flowhelds.  The  results  obtained  by  implicit 
schemes  are  therefore  shown  below. 

The  code  has  both  2D  and  axisymmetric  capabilities,  with  the  latter  one  used  in  this 
work.  A  second  order  axisymmetric  spatial  discretization  is  used  along  with  the  hrst  order 
velocity  {vx,  Vr,  vg)  discretization.  Here  Vx  is  the  velocity  component  along  x  direction,  Vr 
the  radial  component,  and  vg  completes  the  right  handed  system.  Note  that  the  coordinate 
system  used  in  the  hnite  volume  scheme  is  different  from  that  used  in  the  statistical  BGK 
method,  where  Cartesian  components,  Vx,  Vy,  Vz,  were  used  instead.  A  supersonic  inhow 
condition  is  used  at  the  nozzle  throat,  and  vacuum  outhow  conditions  are  set  at  the  outer 
boundaries.  Fully  diffuse  rehection  with  complete  energy  accommodation  is  applied  at  the 
nozzle  surface.  The  spatial  grid  convergence  was  achieved  by  increasing  the  number  of  nodes 
from  3,600  to  over  17,000.  The  convergence  on  the  velocity  grid  was  also  studied,  with  the 
number  of  grid  points  ranging  from  (20,10,18)  to  (30,35,50). 

III.  Results  and  Discussion 

A  gas  how  of  argon  through  a  conical  nozzle  expanding  into  a  vacuum  is  considered 
in  this  work.  The  diverging  portion  of  the  nozzle  is  modeled,  and  its  geometry  is  taken 
from  Ref.  30.  The  nozzle  throat  diameter  is  2.5  mm,  the  length  of  the  diverging  part  is 
50.7  mm,  and  the  half- angle  is  20°.  The  surface  temperature  of  the  nozzle  is  assumed  to  be 
300  K.  The  numerical  results  are  obtained  for  two  throat-diameter  based  Reynolds  numbers 
of  1,230  and  12,300,  with  a  stagnation  temperature  of  333.33  K.  The  respective  how  Knudsen 
numbers  are  1.59x10“^  and  1.59x10““^,  where  the  mean  free  path  has  been  calculated  for  the 
stagnation  conditions.  In  all  numerical  approaches,  the  computational  domain  starts  at  the 
nozzle  throat  and  covers  the  entire  diverging  part  of  the  nozzle,  as  well  as  a  small  part  of 
the  plume  to  avoid  the  inhuence  of  the  downstream  boundary  conditions. Four  diherent 
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numerical  approaches  are  used,  as  explained  in  the  last  section. 

We  compare  the  results  of  the  statistical  BGK/ES-BGK  methods  with  those  of  the  finite 
volume  solution  to  the  same  equation  and  the  DSMG  method.  In  order  to  reduce  the 
statistical  error  so  that  small  differences  among  the  different  kinetic  methods  may  be  studied, 
a  large  number  of  simulated  particles  were  used  such  that  the  minimum  number  of  particles 
per  cell  is  about  150.  This  unusually  large  number  of  particles  per  cell  also  provides  a 
better  statistical  approximation  of  the  local  instantaneous  cell  based  temperature  and  more 
accurate  conservation  of  momentum  and  energy.^® 

Since  the  specific  advantages  of  statistical  BGK  and  ES-BGK  methods  over  other  meth¬ 
ods  is  of  interest,  a  comprehensive  sensitivity  study  is  carried  out  with  respect  to  cell 
size/time  step  and  number  of  particles  in  a  cell. 

III. A.  Low  pressure  Case  I 

Results  for  the  low  pressure  Gase  I  were  obtained  by  the  baseline  DSMG  Smile^^  code,  a 
modihed  version  that  implements  the  statistical  BGK  and  ES-BGK  scheme,  a  hnite  volume 
ES-BGK  solver,  and  a  Navier-Stokes  solver. The  DSMG  method  is  understood  to  provide 
the  most  accurate  solution  for  this  case,  as  the  conditions  are  well  within  the  range  where 
a  mesh  converged  DSMG  solution  is  obtainable  and  all  of  the  strict  DSMG  method  criteria 
(as  described  in  Sec.  HlH)  are  satisfied  during  the  calculation.  Since  temperature  is  the 
most  sensitive  parameter  in  BGK  model  due  to  the  constraint  of  Prandtl  number  equal 
to  unity,  we  show  temperature  contour  plots  for  the  different  methods  for  the  purpose  of 
comparison.  In  all  the  figures,  temperature  is  normalized  by  the  stagnation  temperature  of 
the  flow  (333.33  K),  flow  velocity  is  normalized  by  the  sonic  velocity  at  the  nozzle  throat 
(294.55  m/s),  density  is  normalized  by  the  density  at  the  nozzle  throat  (0.013989  kg/m^),  x 
and  y  coordinates  are  normalized  by  the  throat  diameter  (2.5  mm)  of  the  nozzle,  x/D  =  0 
corresponds  to  the  nozzle  throat  and  y/D  =  {)  corresponds  to  the  nozzle  centerline. 

The  general  flow  physics  of  supersonic  expansions  to  vacuum  are  illustrated  in  Fig.  ^ 
which  shows  a  comparison  of  temperature  contours  predicted  by  the  DSMG,  ES-BGK,  and 
NS  approaches.  All  gas  dynamic  methods  show  the  typical  structure  of  an  inviscid  flow  core 
as  well  as  the  development  of  a  thick,  viscous  boundary  layer.  The  predicted  flow  structure 
for  the  ES-BGK  agrees  best  with  the  DSMG,  whereas,  the  NS  shows  the  poorest  agreement  of 
all  the  gas  dynamic  techniques.  The  other  gas  dynamic  techniques  show  similar  qualitative 
agreement  with  DSMG  compared  to  the  ES-BGK  method,  and,  subtle  differences  will  be 
illustrated  in  linear  prohles  discussed  below. 

To  quantify  the  differences  among  the  different  numerical  methods,  selected  flow  param¬ 
eters  are  presented  along  the  nozzle  centerline  (Fig.|21)  and  across  the  nozzle  at  the  exit  plane 
(Fig.  in}  for  all  of  the  methods.  Figure  El  shows  that  all  gas  dynamic  approaches  agree  well 
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within  the  core  of  the  flow  predicting  virtually  identical  density,  velocity,  and  temperature 
profiles  for  similar  to  (but  not  the  same  as)  isentropic  expansions.  Important  difference 
among  the  techniques,  however,  can  be  seen  in  Fig.  IHlthe  most  likely  portion  of  the  flow 
where  differences  could  occur,  i.e.,  the  expansion  of  the  boundary  layer  over  the  nozzle  lip. 
Figure  El  shows  that  the  statistical  BGK  scheme  deviates  slightly  from  DSMC  close  to  the 
wall  at  the  nozzle  exit  plane.  The  results  obtained  by  the  statistical  ES-BGK  scheme  match 
well  with  those  of  DSMG  and  FV  ES-BGK  method  along  the  nozzle  centerline  and  exit 
plane.  The  FV  ES-BGK  scheme  agrees  well  with  the  DSMG  method,  with  small  deviations 
from  it  close  to  the  nozzle  exit  towards  the  nozzle  wall.  Generally,  FV  ES-BGK  method  can 
be  expected  to  be  close  to  NS  for  larger  Reynolds  number,  and  somewhere  between  DSMG 
and  NS  for  smaller  Reynolds  number.  The  reason  for  this  is  the  limitations  posed  by  the 
ES-BGK  equation  itself,  wherein  the  collision  term  limits  ES-BGK  applicability  to  highly 
nonequilibrium  flows.  The  statistical  solution  of  the  model  kinetic  equation  (BGK/ES-BGK) 
successfully  captures  the  boundary  layer  and  nozzle  core  flow  as  well.  The  NS  solver  pro¬ 
duces  results  accurate  in  the  core  flow  but  strongly  differs  from  the  ES-BGK,  DSMG  and 
particle  BGK  throughout  the  boundary  layer  even  with  the  imposition  of  velocity  slip  and 
temperature  jump  boundary  conditions  at  the  nozzle  wall.  It  can  be  stated  that  NS  pre¬ 
diction  in  the  rarefied  flow  regime  is  inadequate  and  should  not  be  used  when  an  accurate 
solution  of  the  flow  over  the  nozzle  lip  is  required,  such  as  in  contamination  problems.  These 
observations  are  consistent  with  the  modeled  boundary  layer  profiles  obtained  by  the  DSMG, 
the  statistical  BGK/ES-BGK,  and  NS  methods,  as  will  be  discussed  in  Sec.  IIII.Gl 

Having  compared  the  supersonic  flow  expansion  inside  the  nozzle,  we  now  discuss  the  flow 
results  in  the  backflow  region  downstream  of  the  nozzle  exit.  Figure  0] shows  the  variation  of 
axial  velocity  and  temperature  across  a  plane  normal  to  the  nozzle  centerline  at  x/D  of  ~0.85 
downstream  of  the  nozzle  exit  plane.  In  the  same  figures,  the  local  Knudsen  number  based 
on  density  gradient  {Kn  =  ^),  is  also  shown.  It  can  be  clearly  seen  that  the  local  Knudsen 
number  shows  a  large  increase  at  the  nozzle  lip  (y/D=~4),  increasing  to  a  value  greater 
than  unity,  and,  thereby  making  the  NS  equations  invalid  for  modeling  the  back  flow  region. 
From  velocity  and  temperature  profiles,  as  shown  in  Fig.0J  it  can  be  seen  that  the  statistical 
ES-BGK  method  shows  very  good  agreement  with  the  benchmark  DSMG  method  for  both 
axial  velocity  and  temperature  (with  a  maximum  deviation  of  less  than  5%),  while  the  NS 
method  shows  significant  deviation  from  the  DSMG  method.  The  disagreement  between  the 
NS  and  DSMG  method  is  more  than  30%,  thereby  making  the  use  of  NS  method  for  such 
multiscale  flow  applications  questionable.  Using  the  back  fiowfield  data,  we  computed  the 
mass  flow  rate  at  the  exit  plane  of  the  nozzle  in  the  direction  opposite  to  that  of  the  nozzle 
core  flow,  and  found  that  the  statistical  ES-BGK  and  DSMG  methods  provide  a  value  of 
~7.5xl0“®  kg/s,  while  the  NS  solution  does  not  predict  any  mass  flow  rate  in  the  backward 
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direction.  While  it  is  true  that  the  backwards  mass  flow  rate  is  very  small  compared  to  the 
overall  mass  flow  rate  of  20.2x10“®  kg/s  through  the  nozzle,  the  flow  rate  is  such  that  it 
cannot  be  ignored  from  the  standpoint  of  spacecraft  contamination. 

All  the  aforesaid  studies  with  respect  to  the  particle  BGK  and  ES-BGK  methods  were 
carried  out  with  a  cell  size  of  24  pm,  a  time  step  of  2xl0“®s  (satisfying  the  requirements  of 
the  DSMC  method  described  in  Sec.  HlH)  and  about  150  number  of  particles  per  cell  so  as 
to  eliminate  the  statistical  error.  Having  demonstrated  the  effectiveness  of  the  model  kinetic 
equations  in  simulating  a  transitional  micro  nozzle  flow,  we  discuss  the  numerical  efficiency 
of  particle  BGK/ES-BGK  methods  in  comparison  to  DSMC.  To  this  end,  sensitivity  studies 
were  carried  out,  wherein  cell  size,  time  step  and  number  of  particles  in  a  cell  are  varied 
to  reach  an  accurate  and  computationally  most  efficient  solution.  It  was  found  that  the 
mesh  independent  solution  for  particle  BGK  and  ES-BGK  methods  can  be  obtained  with  a 
cell  size  of  80  pm,  a  corresponding  time  step  of  1.5xl0“®s  and  with  an  average  number  of 
about  40  particles  in  a  cell  rather  than  150.  It  took  about  6.5  hours  for  statistical  BGK  and 
ES-BGK  methods  to  reach  a  converged  state  using  20  Intel  3GHz  Xeon  processors.  In  the 
use  of  DSMC,  a  mesh  independent  converged  solution  could  be  obtained  with  a  maximum 
cell  size  of  40  pm,  a  time  step  of  lxl0“®s,  and  with  15  particles  in  a  cell.  The  solution 
required  approximately  10  hours  using  the  same  computer  hardware.  We  note  that  with 
these  numerical  parameters,  the  solutions  obtained  with  the  particle  BGK/ES-BGK  method 
agree  well  with  those  obtained  from  DSMC. 

Table  ^  shows  the  comparison  of  computational  time  required  by  the  various  methods  to 
reach  a  mesh  independent  converged  state  for  the  low  pressure  case.  The  performance  data 
for  DSMC  and  statistical  methods  were  obtained  using  Intel  3GHz  Xeon  processors,  while 
Intel  Itanium  2(1.6  GHz)  processors  were  used  for  the  FV  ES-BGK  and  NS  solver.  It  is 
evident  that  the  particle  BGK  and  ES-BGK  methods  require  less  CPU  time  than  the  DSMC 
method  without  incurring  any  loss  in  the  accuracy  of  the  solution.  It  can  also  be  inferred  that 
the  cell  size  could  be  larger  in  the  statistical  BGK  and  ES-BGK  methods  compared  to  the 
DSMC  method,  which  would  allow  us  to  use  these  methods  for  higher  pressure  cases,  where 
DSMC  cannot  be  used  due  to  the  high  computational  cost.  It  is  interesting  to  note  that 
the  convergence  process  is  quite  different  for  DSMC  and  statistical  BGK/ES-BGK  methods. 
The  DSMC  method  tends  to  reach  steady  state  faster,  but  requires  more  computational 
effort  to  collect  sufficient  information  for  the  solution  to  be  smooth  (i.e.,  sampling).  On  the 
other  hand,  statistical  BGK/ES-BGK  methods  reach  steady  state  slower  (perhaps  due  to 
the  ’’history”  of  the  macroparameter  sampling  procedure  which  defines  the  local  Maxwellian 
distribution  function),  but  smoothness  of  the  results  is  achieved  earlier  than  in  the  case  of 
DSMC  (probably  for  the  same  reason).  The  most  time  efficient  method  is  the  NS  method, 
however,  at  the  cost  of  accuracy  in  the  rarefied  portion  of  the  flow. 
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In  summary  for  this  case,  we  can  conclude  that  the  statistical  and  FV  ES-BGK  results 
are  in  good  agreement  with  the  DSMC  solution  in  the  entire  diverging  part  of  the  nozzle. 
The  statistical  solution  of  the  model  kinetic  equation  accurately  captures  both  the  boundary 
layer  and  the  nozzle  coreflow  with  a  significant  saving  in  computational  time  over  DSMC 
method.  Even  though  the  velocity  slip  and  temperature  jump  boundary  conditions  were 
used  at  the  wall  in  the  case  of  NS,  it  fails  to  model  the  rarefaction  effects  and  predict  the 
proper  macroparameter  values  in  the  boundary  layer  and  at  the  nozzle  lip. 

III.B.  High  pressure  Case  II 

Results  for  the  high  pressure  Case  II  were  obtained  with  the  baseline  DSMC  method,  sta¬ 
tistical  BCK  and  ES-BCK  schemes,  FV  ES-BCK  solver,  and  the  NS  solver.  Although  the 
DSMC  method  is  not  usually  applied  for  such  a  case  due  to  signihcant  computational  re¬ 
quirements,  we  were  able  to  obtain  a  mesh  independent  solution  and  satisfy  all  the  DSMC 
numerical  requirements.  Similar  to  the  low  pressure  case.  Fig.  El  presents  a  comparison  of 
the  predicted  spatial  temperature  distribution  predicted  by  the  DSMC,  ES-BCK,  and  NS 
approaches.  This  and  subsequent  figures  are  normalized  in  the  same  manner  as  discussed 
in  the  previous  section.  The  flow  exhibits  similar  features,  compared  to  the  lower  pressure 
case,  of  an  inviscid  core  and  a  thick  boundary  layer.  The  higher  stagnation  pressure,  how¬ 
ever,  causes  the  boundary  layer  to  be  thinner  compared  to  the  lower  pressure  case  as  can  be 
seen  by  comparison  of  Figs.  Elands  Also,  since  the  present  case  is  for  a  pressure  10  times 
higher  than  for  the  low  pressure  case,  the  agreement  between  the  solutions  of  the  NS  and 
DSMC  method  is  better  than  in  the  lower  pressure  case,  as  expected.  However  as  before, 
the  agreement  slightly  worsens  in  the  relatively  rarehed  portions  of  the  flow,  near  the  nozzle 
exit  and  close  to  the  boundary  layer. 

To  provide  a  quantitative  evaluation,  we  now  compare  results  of  all  of  the  numerical 
schemes  in  several  linear  plots  presenting  the  velocity  and  temperature  profiles  along  the 
nozzle  center  line  (Fig.  E}  and  at  the  exit  plane  of  the  nozzle  (Fig.  [7j).  Figures  El  and  [7|  also 
present  results  of  the  eDSMC^*^  solution  for  this  high  pressure  case.  The  eDSMC  calculations 
were  computed  using  0.85  million  cells  with  5  million  particles.  Figure  El  shows  that  the 
DSMC,  NS,  FV  ES-BCK,  statistical  BCK,  ES-BCK,  and  eDSMC  predicted  profiles  for  both 
velocity  and  temperature  are  close  along  the  nozzle  centerline.  Because  the  nozzle  throat 
region  is  modeled  in  our  studies  as  a  sharp  corner,  the  flow  develops  a  weak  shock  wave  which 
is  then  reflected  from  the  centerline.  It  can  be  seen  that  the  DSMC,  FV  ES-BCK  and  NS 
methods  were  able  to  predict  this  feature  of  the  flow,  but  the  statistical  BCK  and  ES-BCK 
methods  slightly  overpredict  viscosity,  diffusing  this  feature  through  the  weak  shock  wave. 
The  lack  of  prediction  of  the  weak  shock  wave  also  accounts  for  some  of  the  difference  in  the 
inviscid  core  seen  between  the  DSMC  and  ES-BCK  contours  shown  in  Fig.  El  Across  the 
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nozzle  exit  plane,  as  in  case  I,  there  is  more  variation  in  the  velocity  and  temperature  prohles 
among  the  different  methods,  as  seen  in  Fig.  [7|  The  hgure  shows  that  the  eDSMC  method 
predicts  a  velocity  and  thermal  boundary  layer  profile  signihcantly  different  from  the  other 
methods,  with  even  small  differences  observed  in  the  core  region  {y/D  <  2)  temperature 
prohle.  This  disagreement  is  consistent  with  that  observed  in  earlier  work  even  though  this 
high  pressure  case  is  a  factor  of  ten  more  dense  than  was  considered  in  Ref.  20.  The  results 
across  the  nozzle  exit  plane  show  that  there  is  still  some  difference  between  the  NS  and 
DSMC  predictions  in  both  the  temperature  and  flow  velocity.  As  in  the  lower  pressure  case, 
this  difference  is  explained  by  the  limitations  of  the  velocity  slip  and  temperature  jump 
boundary  condition  used  in  the  NS.  The  remaining  numerical  methods  all  agree  well  with 
DSMC. 

In  a  manner  similar  to  the  low  pressure  case,  we  now  present  the  flowfield  data  in  the 
back  flow  region  for  Case  II.  Figure IHlshow  the  variation  of  the  x  component  of  velocity  and 
temperature  across  a  plane  normal  to  the  nozzle  centerline  at  x/D  =~  0.85  downstream  of 
the  nozzle  exit  plane.  Again,  the  variation  of  local  Knudsen  number  based  on  the  density 
gradient,  is  shown  in  the  hgure  and  it  can  be  seen  that  the  variation  is  large  in  the  nozzle 
lip  region,  even  for  a  high  pressure  case.  Similar  to  the  low  pressure  case,  the  solution 
predicted  by  the  NS  equations  deviates  from  that  of  DSMC.  The  hgure  shows  that  although 
there  are  regions  with  Knudsen  number  well  within  the  continuum,  deviation  between  the 
two  solutions  is  signihcant,  particularly,  in  the  region  where  Kn  <  0.1.  The  disagreement 
between  the  NS  and  DSMC  is  as  much  as  20%  in  both  velocity  and  temperature,  whereas  the 
statistical  ES-BGK  method  shows  only  a  discrepancy  of  less  than  4%  from  the  DSMC,  which 
is  within  the  numerical  accuracy.  Once  again,  the  back  flow  rate,  calculated  at  the  nozzle 
exit  plane,  predicted  by  the  NS  solution  is  zero,  whereas  both  DSMC  and  ES-BGK  methods 
predict  a  value  of  ~  18x10“®  kg/s  (compared  to  the  total  mass  flow  rate  of  202.2x10“®  kg/s 
through  the  nozzle). 

Similar  to  the  low  pressure  case,  we  carried  out  a  sensitivity  study  for  this  case  to  demon¬ 
strate  the  numerical  efficiency  of  the  methods.  Again,  the  performance  data  for  the  DSMC 
and  statistical  BGK/ES-BGK  methods  was  obtained  using  Intel  3GHz  Xeon  processors, 
while  Intel  Itanium  2(1.6  GHz)  processors  were  used  for  the  FV  ES-BGK  and  NS  solver. 
It  was  found  that  it  takes  less  time  to  reach  a  mesh  independent  converged  solution  using 
the  particle  BGK  and  ES-BGK  methods  compared  to  a  DSMC  converged  solution.  Table  |21 
shows  the  comparison  of  computational  time  required  by  the  various  methods  for  the  high 
pressure  case.  It  should  be  noted  that  the  particle  BGK  and  ES-BGK  methods  take  less 
time  compared  to  DSMC  because  the  cell  size  requirement  is  relaxed  in  the  BGK  methods, 
although  these  methods  require  more  particles  per  cell  than  DSMC.  This  case  is  perhaps 
the  limiting  maximum  pressure  case  solvable  by  the  DSMC  method,  with  the  present  com- 
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putational  power  making  it  unfeasible  to  extend  to  higher  pressures.  Since  the  cell  size 
requirement  can  be  relaxed  for  the  statistical  BGK  and  ES-BGK  methods,  they  can  be  used 
for  higher  pressure  cases.  The  FV  ES-BGK  method  takes  a  large  amount  of  time,  although, 
its  results  agree  well  with  the  DSMC  method.  The  eDSMG  method,  although  fast  and  com¬ 
parable  with  the  NS  solver,  does  not  capture  the  viscosity  effects  in  the  boundary  layer  as 
was  suggested  by  earlier  results  discussed  in  Ref.  20.  However  the  solution  in  the  inviscid 
core  of  the  flow  is  remarkably  close  to  the  solutions  obtained  by  the  other  methods.  The  NS 
solver  again  turns  out  to  be  the  most  time  efficient  method,  however  at  the  cost  of  accuracy 
in  the  rarehed  portions  of  the  boundary  layer. 

III.C.  Comparison  of  Averaged  Flow  Quantities  for  Low  and  High  Pressure 
Cases 

We  now  discuss  the  averaged  flow  quantities,  such  as  boundary  layer  thickness,  displacement 
thickness  and  thrust,  derived  from  the  flowheld  data  obtained  by  the  different  methods. 
The  boundary  layer  thickness,  a  function  of  distance  along  the  centerline,  is  dehned  as 
that  distance  from  the  nozzle  wall  where  the  local  flow  velocity  reaches  0.99  of  the  axial 
velocity.^^  Note  that  the  presence  of  rarefaction  causes  the  boundary  layer  thickness  to 
deviate  from  the  traditional  slender  body  continuum  square  root  of  axial  distance  behavior. 
Figure  El  compares  the  boundary  layer  prohles  obtained  by  the  DSMC,  NS  and  the  statistical 
BGK  and  ES-BGK  methods  for  the  low  and  high  pressure  cases.  It  can  be  seen  that  the 
agreement  between  the  statistical  BGK/ES-BGK  and  DSMC  methods  is  good,  while,  there 
is  a  noticeable  difference  between  the  NS  and  DSMC  values.  Table Olshows  a  comparison  of 
the  boundary  layer  thickness  at  the  nozzle  exit  obtained  by  the  aforesaid  methods  for  the 
low  and  high  pressure  cases.  The  percentage  difference  between  the  values  of  boundary  layer 
thickness  (at  the  nozzle  exit)  obtained  from  the  NS  and  DSMC  methods  is  on  average  about 
8%  for  both  pressure  cases.  The  statistical  BGK  and  ES-BGK  methods  show  considerably 
good  agreement  DSMC  with  a  maximum  deviation  of  about  4.8%  for  the  low  and  high 
pressure  cases. 

The  displacement  thickness  is  dehned  as  the  distance  the  undisturbed  core  how  is  dis¬ 
placed  from  the  boundary  by  a  stagnant  layer  that  removes  the  same  mass  how  from  the 
how  held  as  the  actual  boundary  layer,  i.e.,  it  is  the  thickness  of  a  zero  velocity  layer  that 
has  the  same  mass-how  defect  as  the  actual  boundary  layer. Table|3]shows  a  comparison  of 
the  displacement  thickness  at  the  nozzle  exit  for  the  aforesaid  methods  for  the  low  and  high 
pressure  cases.  It  is  worth  noticing  that  the  NS  solution  results  in  a  signihcantly  diherent 
displacement  thickness  values  relative  to  DSMC  as  compared  to  the  other  methods.  The 
percentage  diherence  between  the  NS  and  DSMC  values  is  greater  than  ~  40%  for  both  low 
and  high  pressure  cases.  On  the  other  hand,  the  statistical  BGK/ES-BGK  methods  show 
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good  agreement  with  the  DSMC  methods  with  a  maximum  disagreement  of  ~  4.5%  for  the 
low  and  high  pressure  cases.  From  Tables  El  and  EJ  it  can  be  inferred  that  the  statistical 
methods  resolve  the  boundary  layer  quite  well,  which  in  turn  would  suggest  that  they  will 
be  useful  for  accurately  predicting  other  derived  flow  parameters  such  as  thrust  as  discussed 
below. 

Using  the  flow  held  data,  thrust  was  computed  for  the  DSMC,  statistical  BGK/ES-BGK 
and  NS  methods.  Table  El  compares  the  thrust  values  obtained  for  the  different  methods  using 
their  respective  how  held  data  for  both  low  and  high  pressure  cases.  It  should  be  noted  that 
the  percentage  diherence  between  the  thrust  values  obtained  from  the  diherent  methods  is 
less  than  3%  for  both  low  pressure  and  high  pressure  cases.  The  maximum  disagreement  of 
3%  is  observed  between  NS  and  DSMC  solutions  for  the  low  pressure  case.  The  thrust  is  a 
spatially  averaged  quantity  that  includes  the  core  and  boundary  layer  portions  of  the  how. 
Since  the  inviscid  core  how,  particularly  at  higher  Reynolds  numbers,  spatially  dominates  the 
average,  the  thrust  agreement  between  diherent  methods  is  better  than  that  computed  for 
the  boundary  layer  thickness.  The  statistical  methods  provide  an  accurate  alternative  to  NS 
for  predicting  micro- nozzle  thrust,  but,  at  a  signihcant  computational  cost.  However,  even 
though  the  NS  solution  provides  a  reasonably  accurate  prediction  of  thrust,  it  fails  to  model 
the  back  how  correctly  as  was  shown  in  Figs.  E]  and  El  for  the  low  and  high  pressure  cases, 
respectively.  The  numerical  advantage  associated  with  the  statistical  ES-BGK  methods, 
compared  to  DSMG,  together  with  its  accuracy,  makes  it  a  preferable  method  for  modeling 
multi-scale  hows,  involving  contaminant  back  how,  where  in  particular,  DSMG  method  may 
not  be  usable  at  high  pressures. 


IV.  Conclusions 

Argon  how  through  a  conical  nozzle  was  studied  for  two  Reynolds  numbers  of  1,230  and 
12,300,  using  four  diherent  approaches.  These  include  a  continuum  approach,  solution  of 
Navier-Stokes  equations,  and  three  kinetic  approaches,  the  DSMG  method,  and  statistical 
and  deterministic  methods  for  the  BGK/ES-BGK  equations.  Analyses  of  the  accuracy  of  the 
approaches  and  their  numerical  efficiency  were  conducted  for  supersonic  expansions  to  vac¬ 
uum  nozzle  hows.  These  cases  were  selected  because  they  span  the  range  of  conditions  found 
in  micro-nozzle  hows.  Since  the  DSMG  numerical  parameters  were  chosen  to  ensure  that  its 
solutions  were  independent  of  the  numerical  parameters  for  both  the  high  and  low  pressure 
cases,  it  was  considered  truth  in  the  comparison  with  the  other  gas  dynamic  approaches. 
Several  conclusions  can  be  drawn  from  the  results  of  the  computations.  The  statistical  and 
hnite  volume  solution  of  the  BGK  and  ES-BGK  equations  are  in  a  good  agreement  with  the 
DSMG  method  in  the  entire  computational  domain  for  both  Reynolds  numbers.  Statistical 


Distribution  A:  Approved  for  public  release;  distribution  unlimited 


BGK  and  ES-BGK  methods  save  computational  time  compared  to  DSMG  without  incurring 
any  loss  in  the  accuracy  of  the  solution.  Also  it  was  found  that  the  statistical  BGK/ES- 
BGK  methods  did  not  require  the  use  of  the  strict  DSMG  numerical  criteria,  as  described  in 
Sec.lThAI  and  therefore  are  more  likely  candidates  for  use  in  multi-scale,  high  pressure  cases, 
where  the  use  of  the  DSMG  method  is  either  not  possible  or  not  practical.  The  eDSMG 
method  was  not  able  to  reproduce  the  viscosity  effects  in  the  boundary  layer,  although, 
because  the  solution  in  the  inviscid  core  is  good  it  could  be  used  in  a  hybrid  approach. 

The  Navier-Stokes  solutions  are  computationally  the  most  tractable  and  are  in  a  good 
agreement  with  the  DSMG  results  in  the  higher  density  core  portion  of  the  flow  where 
rarefaction  effects  are  small.  In  the  boundary  layer,  however,  even  though  the  velocity 
slip  and  temperature  jump  boundary  conditions  were  used,  there  is  a  noticeable  difference 
between  the  NS  and  DSMG  solutions.  Although  the  NS  method  results  in  a  reasonably 
accurate  thrust  prediction,  it  fails  to  model  correctly  the  back  flow  beyond  the  nozzle  exit, 
which  is  important  when  contamination  of  mission  critical  systems  from  onboard  spacecraft 
nozzle  systems  must  be  avoided. 
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Table  1.  Efficiency  of  the  methods,  Low  Pressure  Case  I 


Method 

DSMG 

Part  BGK 

Part  ES-BGK 

FV  ES-BGK 

NS 

Time  (GPUH) 

200 

130 

130 

400 

<1 

Num.  of  Particles  (million) 

15 

4.67 

4.67 

- 

- 

Num.  Gells 

3,000,000 

120,000 

120,000 

3,600 

3,600 

Table  2.  Efficiency  of  the  methods,  High  Pressure  Case  II 


Method 

DSMG 

Part  BGK 

Part  ES-BGK 

FV  ES-BGK“ 

NS 

Time  (GPUH)^ 

1000 

840 

840 

5000 

<3 

Number  of  Particles  (million) 

75 

50 

50 

- 

- 

Number  of  Gells 

15,000,000 

1,260,000 

1,260,000 

17,00 

14,400 

“FV  ES-BGK  times  are  for  spatial-velocity  grid  converged  results 

^Tlie  times  reported  for  the  first  three  methods  are  for  a  3GHz  Xeon  processor, 

while  the  last  two  are  for  a  Intel  Itanium  2(1.6  GHz)  processor. 


Table  3.  Comparison  of  boundary  layer  thickness  (mm)  at  the  nozzle  exit,  Case  I,  II 


Gase 

DSMG 

Particle  BGK 

Particle  ES-BGK 

NS 

Re=1230  (Gase  I) 

7.73 

8.09 

7.83 

7.20 

Re=12300  (Gase  II) 

5.96 

6.25 

6.14 

5.40 

Table  4.  Comparison  of  displacement  thickness  (mm)  at  the  nozzle  exit,  Case  I,  II 


Gase 

DSMG 

Particle  BGK 

Particle  ES-BGK 

NS 

Re=1230  (Gase  I) 

6.47 

6.19 

6.42 

3.66 

Re=12300  (Gase  II) 

3.34 

3.18 

3.23 

1.87 
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Table  5.  Comparison  of  nozzle  thrust  (mN),  Case  I,  II 


Case 

DSMC 

Particle  BGK 

Particle  ES-BGK 

NS 

Re=1230  (Case  I) 

9.85 

9.77 

9.81 

9.55 

Re=12300  (Case  II) 

1.12 

1.10 

1.11 

1.10 

(a)  Statistical  ES-BGK  vs.  DSMC  (b)  NS  vs.  DSMC 

Figure  1.  Case  I,  comparison  of  temperature  [K]  contours  for  different  flow  methods. 


(a)  Density  (b)  X  component  velocity  (c)  temperatnre  (K) 

Figure  2.  Case  I,  comparison  of  nozzle  centerline  macro- parameters  for  different  gas  dynamic  techniques. 
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(a)  X  component  velocity  (b)  Temperature 

Figure  3.  Case  I,  comparison  of  velocity  and  temperature  profiles  across  the  nozzle  exit  for  different  gas 
dynamic  approaches. 


Knudsen  Number 

10''  10^  10''  10'  10" 


Knudsen  Number 

10''  10"  10"  10'  10" 


(a)  X  component  velocity  (b)  Temperature 

Figure  4.  Case  I,  comparison  of  velocity,  temperature,  and  local  Knudsen  number  profiles  across  a  plane 
normal  to  the  nozzle  centerline  at  a  location  of  x/D~0.85  beyond  the  nozzle  exit  for  different  gas  dynamic 
approaches. 


Distribution  A:  Approved  for  public  release;  distribution  unlimited 


n/n 


x/D  x/D 

(a)  Statistical  ES-BGK  vs.DSMC  (b)  NS  vs.  DSMC 

Figure  5.  Case  II,  comparison  of  temperature  [K]  contours  for  different  flow  methods. 


(a)  X  component  of  velocity 


(b)  Temperature 


Figure  6.  Case  II,  comparison  of  nozzle  centerline  macro-parameters  for  different  gas  dynamic  techniques. 
Case  II,  X  component  velocity  profile  along  the  nozzle  centerline 
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(a)  X  component  of  velocity  (b)  Temperature 

Figure  7.  Case  II,  comparison  of  velocity  and  temperature  profiles  across  the  nozzle  exit  for  different  gas 
dynamic  approaches. 


Knudsen  Number 


Knudsen  Number 

lO"*  10^  10^  10'  10“ 


T/To 


(a)  X  component  velocity  (b)  Temperature 

Figure  8.  Case  II,  comparison  of  velocity,  temperature,  and  local  Knudsen  number  profiles  across  a  plane 
normal  to  the  nozzle  centerline  at  a  location  of  x/D^O.85  beyond  the  nozzle  exit. 
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(a)  Case  1  (b)  Case  II 

Figure  9.  Comparison  of  growth  of  the  boundary  layer  thickness  along  the  nozzle  for  the  low  and  high  pressure 
cases. 
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